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ABSTRACT 

We carried out three-dimensional global resistive magnetohydrodynamic 
(MHD) simulations of the cooling instability in optically thin hot black hole ac- 
cretion flows by assuming bremsstrahlung cooling. General relativistic effects are 
simulated by using the pseudo-Newtonian potential. Cooling instability grows 
when the density of the accretion disk becomes sufficiently large. We found 
that as the instability grows the accretion flow changes from an optically thin, 
hot, gas pressure-supported state (low/hard state) to a cooler, magnetically sup- 
ported, quasi-steady state. During this transition, magnetic pressure exceeds the 
gas pressure because the disk shrinks in the vertical direction almost conserving 
the toroidal magnetic flux. Since further vertical contraction of the disk is sup- 
pressed by magnetic pressure, the cool disk stays in an optically thin, spectrally 
hard state. In the magnetically supported disk, the heating rate balances with 
the radiative cooling rate. The magnetically supported disk exists for time scale 
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much longer than the thermal time scale and comparable to the accretion time 
scale. 

We examined the stability of the magnetically supported disk analytically, 
assuming that the toroidal magnetic flux is conserved, and found it thermally 
and secularly stable. Our findings may explain why black hole candidates stay 
in luminous, X-ray hard state even when their luminosity exceeds the threshold 
for the onset of the cooling instability. 

Subject headings: accretion, accretion disks — black hole physics — magnetohydro- 
dynamics:MHD 

1. Introduction 

State transitions between a low/hard state and a high/soft state are observed in Galactic 
black hole candidates such as Cyg X-l (e.g., Holt et al. (1975); Oda (1977); Zdziarski et al. 
(2002)), GS 1124-68 (e.g., Kitamoto et al. (1992); Ebisawa et al. (1994)), GRS 1915+105 
(e.g., Belloni et al. (2000); Ueda et al. (2002)) and GX 339-4 (e.g., Miyamoto et al. (1991), 
1995; Zdziarski et al. (2004); Homan et al. (2005)). Recent RXTE observations of state 
transitions in binary black hole candidates and their relation to jet ejections are reviewed by 
Fender, Belloni & Gallo (2004) and Remillard (2005). The low/hard state is characterized by 
violent X-ray fluctuations and absence of the soft blackbody component in its spectrum. In 
this state, mass accretes to the black hole as an optically thin, advection dominated accretion 
flow (e.g., Ichimaru (1977); Narayan & Yi (1994), 1995). The spectrum of a high/soft state is 
dominated by blackbody-like component emitted from an optically thick, geometrically thin 
disk (e.g., Shakura & Sunyaev (1973)). Recent X-ray observations of outbursts of GX 339-4 
indicate that the source stays in the spectrally hard and luminous intermediate state during 
the transition from low/hard state to the high/soft state (Homan et al. (2005); Remillard 
(2005)). It has been suggested that relativistic ballistic ejections occur during the transition 
from this intermediate state to the optically thick, soft state (Fender et al. 2004). 

Outbursts of black hole candidates are triggered by the increase of the mass accretion 
rate from the outer radius (Mineshige & Wheeler 1989). In the early stage of an outburst, 
since a hot, optically thin, hard X-ray emitting disk is formed around the black hole, the 
black hole candidate is observed as staying in the low/hard state. The X-ray luminosity of 
the disk increases as the accretion rate increases. Abramowicz et al. (1995) obtained thermal 
equilibrium curves of optically thin accretion disks including both the bremsstrahlung cooling 
and advective cooling (see also Kato et al. (1998)). They showed that when the mass 
accretion rate (and, hence, surface density) exceeds some limit, optically thin disks cannot 
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be in thermal equilibrium because the radiative cooling overwhelms the heating. When this 
limit is reached, therefore, a thermal instability is triggered so that the disk should cool and 
shrink in the vertical direction (Mineshige (1996); Lasota, Narayan, & Yi 1996). 

During this vertical contraction of the disk, magnetic fields in the turbulent disk will 
be amplified due to flux conservation or will be dissipated by magnetic reconnection. In 
the former case, strongly magnetized disk will be formed. The possible existence of such 
low-/? (/? = -P g as/-Pmag < 1) disk was pointed out by Shibata, Tajima & Matsumoto (1990). 
Although magnetic flux can buoyantly escape from the disk by the Parker instability (Parker 
1966), once low-/? disk is formed, the disk stays in the low-/? state because magnetic tension 
prevents the growth of the Parker instability. They called such disks as magnetic cataclysmic 
disks because explosive events will take place when the magnetic energy stored in the disk is 
released. Mineshige, Kusunose & Matsumoto (1995) suggested that low-/? disks may corre- 
spond to the low/hard state. Recently Pariev, Blackman and Boldyrev (2003) constructed 
a model of an optically thick, low-/? disk based on a-prescription of viscosity. However, they 
did not show how the low-/? disks are created. Three-dimensional resistive MHD simulations 
can answer this question. 

Global 3D MHD simulations of radiatively inefficient accretion disks have been carried 
out by several authors (e.g., Matsumoto (1999); Hawley (2000),2001; Machida et al. (2000); 
Armitage & Reynolds (2003); Machida & Matsumoto (2003); Igmenshchev et al. (2003); Kato 
et al. (2004)). It has been shown that the MHD turbulence driven by magneto-rotational 
instability (MRI) enhances the angular momentum transport rate and enables mass accre- 
tion. Initially weak magnetic fields are amplified due to the MRI, and saturates when the 
plasma-/? is about 10. Machida, Nakamura & Matsumoto (2004) showed that the radial 
structure of its innermost region (r < 20r g , where r g is the Schwarzschild radius) agrees 
well with that of one- dimensional steady transonic solutions, in which the radial advection 
and the radial dependence of a are both taken into account. As the accretion rate M in- 
creases, the surface density of the disk S increases. Numerical solutions follow the evolution 
track E oc M expected from conventional theory of optically thin, advection dominated hot 
accretion flows (e.g., Abramowicz et al. (1995)). 

The purpose of this paper is to study the growth of the thermal instability which takes 
place when the density of the optically thin disk becomes sufficiently large. We first simulated 
the formation of an optically thin, turbulent disk without including the radiative cooling. 
Subsequently, we simulated the growth of the thermal instability by including the radiative 
cooling term. 

In section 2, we describe basic equations, initial conditions, and numerical method. 
The results of simulations are given in section 3. Section 4 is devoted for discussion and 
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conclusion. 



2. Numerical Model 

2.1. Basic Equations 

We solved the following resistive MHD equations in a cylindrical coordinate system 
(w,<p,z); 



P 



dt 

dv ^ 



()P + V(pv) = 0, (1) 



-VP - P V(j) + , (2) 
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— = Vx(ux£ 77j) , (3) 

at c 
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where p, P, 0, v, j = cV x B/An, rj, T and 5 are the density, pressure, gravitational po- 
tential, velocity, magnetic field, current density, resistivity, temperature and specific entropy, 
respectively. The specific entropy is expressed as S = C v ln(P/p 7 ) where C v is the specific 
heat capacity and 7 is the specific heat ratio. The Joule heating term and the radiative 
cooling term Q~ ad are included in the energy equation. We assume the anomalous resistivity 
77 = ?7o[max(f d /t> c — 1,0)] 2 (Yokoyama & Shibata 1998) where = j/p is the electron-ion 
drift speed and v c is the threshold above which anomalous resistivity sets in. We assume 
bremsstrahlung cooling Q~ ad = 6.2 x 10 20 p 2 T 1 / 2 erg s _1 cm~ 3 . 

General relativistic effects are simulated by using the pseudo-Newtonian potential <fi = 
—GM/(r — r g ) (Paczyhski & Wiita 1980), where G is the gravitational constant, M is the 
mass of the black hole, and r = (w 2 + z 2 ) 1 / 2 . We neglect the self gravity of the disk. 



2.2. Numerical Methods and Boundary Conditions 

We solved the resistive MHD equations in a cylindrical coordinate system (w, tp, z) by 
using a modified Lax-Wendroff scheme (Rubin & Burstein 1967) with an artificial viscosity 
(Richtmyer & Morton 1967). 
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For normalization, we use units listed in table 1. The units of length and velocity are 
Schwarzschild radius r g and the light speed c. For unit density, we adopt p = l/{K es r g ) = 
8.3 x 1O _7 (M/1OM ) g cm -3 , where K es = 0.4 cm 2 /g is the electron scattering opacity. 
The unit time is to — r g/ c = 10~ 4 (M/10M Q ) sec. The unit temperature is given by T = 
rripC 2 /k-Q = 1.1 x 10 13 K, where m p is the proton mass and /cb is the Boltzmann constant. 
The accretion rate is normalized by Msdd = -^Edd/c 2 , where L-^dd = AncGM/ K es = 1.25 x 
1O 39 (M/1OM ) erg s" 1 is the Eddington luminosity. 

In the energy equation, we use the normalized cooling rate 



Here, T is ion temperature and Q' h is a parameter depending on the electron-ion coupling 
and inverse Compton amplification. In this paper, we assume single temperature plasma 
and adopt Q' h = 1.9 x 10" 4 . 

We solved the energy equation transformed to the conservation form; 



d (\ , , 2 B 



'2 



1 , ,2 lP \ , 1 



-pv" + -L—jv+-E-xB\=-pvV<t>- Q rad 

(6) 

where E = —v x B + Anrjj and ' denotes the normalized quantities. In the following, we 
omit ' from the symbols for simplicity except the radiative cooling term. In this formulation 
of the energy equation, the numerically dissipated kinetic energy and magnetic energy are 
captured as the thermal energy (Hirose et al. 2005). Shock heating is also incorporated 
self-consistently. When magnetic turbulence develops in the disk, the disk plasma is heated 
more efficiently than expected from the Joule heating. However, in the early stage of the 
simulation when the disk is still laminar, since we did not include the explicit heating term 
except the Joule heating, the disk heating is insufficient to balance with the cooling term. 
Thus, we include the radiative cooling term after a quasi-steady turbulent accretion disk is 
formed through accretion from the initial torus. 

The number of grids is (iV ro , iV^, N z ) = (250, 64, 191). The grid size is Azu = Az = 0.1 
for < w, z < 10, and otherwise increases with w and z. The grid size in azimuthal 
direction is A<p = 27r/63. We simulated only the upper half space (z > 0) by applying 
a symmetric boundary condition at the equatorial plane [z — 0). The outer boundaries at 
zu = 132 and at z = 69 are free boundaries where waves can be transmitted. We included the 
full circle (0 < if < 2n) in the simulation region, and applied periodic boundary conditions 
in the azimuthal direction. An absorbing boundary condition is imposed at r = r in = 2 by 
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introducing a damping parameter; 

a .0,(,0- t a„ h ^^). , 7 ) 

Physical quantities q = (p,v, B, P) inside r = r in are re-evaluated by 

q new =q~a(q- %) , (8) 

which means that the deviation from initial values go is artificially reduced with damping 
rate a. Waves propagating inside r = r in are absorbed in the transition region (r in — 5Atu < 
r < r in ). 



2.3. Initial Model 

The initial state of our simulation is an equilibrium torus threaded by weak toroidal 
magnetic field (Okada et al. 1989). Initially the torus is assumed to have a constant specific 
angular momentum L. 

According to Okada et al. (1989), by using the polytropic relation P = Kp 1 at the 
initial state and by assuming 

where /3b is the initial plasma (3 at the initial pressure maximum of the torus (w, z) = (rob, 0), 
and B v is the azimuthal magnetic field, we integrated the equation of motion into a potential 
form; 

*(o7, z) = + ^2 + ^77^ + 2 ( 7 1 a = ^ = constant , (10) 

where t> s = (^P/p) 1 / 2 is the sound speed, t>A = [-B^,/(47rp)] 1//2 is the Alfven speed, and 
\l/b = ^(c^bjO). At cc7 = rob, the rotation speed of the torus L/zuh equals the Keplerian 
velocity. By using equation (10), we obtain the density distribution as 

( max{M/ b -0-i/7(2u; 2 ),O} ( . 

9 'Vb/iT-iiiii+^^ve 1 ']/ 

where pb is the density at (w,z) = (wb,0). Outside the torus, we assumed a hot, isother- 
mal (T = T halo ) spherical halo. The density distribution of the halo is given by p h = 
Phaio ex p[ — (0 — 0t>)/(&B^haio)]) where 0b is the gravitational potential at (va,z) = (zu h ,0). 
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We assume that the radiative cooling is negligible in the coronal region where the density is 
lower than p crit = 10~ 4 pb- 

In this paper, we adopted model parameters zu h = 50, K = P h / = 5 x 1CT 4 , (3 h = 100, 
7 = 5/3, L = (zu h /2) 1 / 2 zu h /(zu h — 1), p ha io = 10~ 4 Pb, Vo = 5 x 10 -4 , and v c = 0.9. The inner 
and outer radii of the initial torus are zu in = 33, and zu ont = 85, respectively. 

3. Numerical Results 

3.1. Evolution of a Thermally Unstable Accretion Disk 

Here we present the results of the numerical simulation for a model with pb = 0.29. 
Figure 1 shows the time evolution of azimuthally averaged density log(p/pb), temperature 
and azimuthal component of magnetic fields. The top panels show the initial state. Figure 
Id, le and If show the stage just before the cooling term is switched on. The initial torus 
deforms itself into an accretion disk by transporting angular momentum through Maxwell 
stress. Azimuthal magnetic fields are not completely random but show sector structures. 

We switched on the radiative cooling term at t = 24100 when the disk becomes quasi- 
steady. At this stage, the Thomson optical thickness of the disk is r d = /t cs S ~ 3. Here, 
the surface density £ is computed by vertically integrating the azimuthally averaged density 
from z = —H d to z = H d) where H d is the half thickness of the disk, which we take 
H d = 0.3^7 + 13.5. The disk is thick for Thomson scattering (r d ~ 3) but effectively optically 
thin (see e.g., Kato et al. (1998)). This state corresponds to the low/hard state observed in 
black hole candidates. 

The third panels from top in figure 1 show the stage after the onset of the cooling 
instability. Geometrically thin, cool, dense disk is created in 15 < zu < 70. Magnetic fields 
in the equatorial region of the vertically shrinking cool disk become almost toroidal and 
well-ordered. In the inner region (zu < 15), the disk still stays in the hot state. 

For comparison, in the bottom panels in figure 1, we show the numerical results at 
t = 32000 for a model without including the radiative cooling term. The disk continues to 
stay in the hot, magnetically turbulent state. Between t = 24000 and t = 32000, the disk 
structure does not change significantly. 

We also confirmed that even when we included the cooling term, if the initial density 
of the torus is small enough (p b < 0.05), the whole disk stays in optically thin hot state 
throughout the simulation (0 < t < 32000). 
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Figure 2a shows the radial distribution of the temperature averaged in the azimuthal 
direction and in the vertical direction (0 < z < 1). Curves show the temperature at t = 
24000, 26000, 30000, and 34000 from top to bottom. This figure indicates that cooling front 
propagates inward. The temperature in the cool region decreases down to 10~ 5 (~ 10 8 K in 
physical units). The region w < 15 stays in the hot state throughout the simulation. The 
inner region stays in the hot state because the accretion rate is not large enough to trigger 
the thermal instability in this region. Note that in conventional models of optically thin, hot 
accretion flows including bremsstrahlung cooling, the critical accretion rate for the onset of 
the cooling instability increases with decreasing the radius (Figure 2 in Abramowicz et al. 
(1995)). 

Figure 2b shows the time evolution of azimuthally averaged surface density, equatorial 
density, equatorial temperature and equatorial plasma (3 at w = 35. Here, surface density 
is integrated in \z\ < 24. Numerical results indicate that after the inclusion of radiative 
cooling, the disk evolves quasi-steadily until t = 27000. Subsequently, the disk shrinks in 
the vertical direction due to the cooling instability. As the cooling instability grows, plasma 
(3 decreases and the disk becomes dominated by magnetic pressure (i.e. (3 < 1). 

Let us compare the cooling time scale t coo \ = -P/Q~ d = (T/To) 1 / 2 /[Q' h (p/ Po)] with other 
time scales. Before the onset of the cooling instability, since T ~ 10 -3 [~ 10 10 K in physical 
units] and p ~ 0.1 [~ 8.3 x 10 _8 (M/10M Q ) g cm -3 ], the cooling time scale is t coo i ~ 1-6 x 10 3 
[~ O.16(M/1OM ) sec in physical units] at w = 35. This time scale is comparable to the 
rotation time scale ~ 27rtd y n ~ 1-8 x 10 3 where tdyn = 1/fi ~ 300 is the dynamical time scale 
at vj = 35. Here fl is the rotation angular speed. The accretion time scale t acc = zu/v m ~ 10 4 
is an order of magnitude larger than the thermal time scale. Numerical results indicate that 
the temperature of the disk decreases from T ~ 10~ 3 to T ~ 10~ 5 in t ~ 6 x 10 3 . This time 
scale (~ 4t cool ) is consistent with the time scale for the growth of the thermal instability. 

Figure 3 shows the isosurface of the plasma (3 at t — 24000 and t = 32000. The blue, 
green, and yellow surfaces show (3 = 100, 10, and 1, respectively. Before the cooling term is 
switched on, \ow-(3 regions ((3 < 1) occupy only a small fraction of the disk. However, after 
the onset of the transition, the low-/? regions fill the disk. The filling factor of the low-/? 
region in \z\ < 1 increases from 0.1 to 0.7 during the transition. 

3.2. Structure of Global Magnetic Fields 

Top panels of figure 4 show the equatorial logarithmic density (color) and magnetic field 
lines projected onto the equatorial plane (gray curves). Before the cooling instability grows 



- 9- 



(figure 4a), magnetic field lines show loosely wound global spiral structure superposed by 
turbulent components. As the disk cools, magnetic fields turn into a tightly wound spiral 
(figure 4b). The equatorial density increases in the outer region (zu > 15), where cooling 
instability grows. Figure 4c shows the equatorial density and magnetic field lines at t — 32000 
for a model without radiative cooling. In this model, the equatorial density remains small 
and the magnetic fields show large amplitude fluctuations. 

In low-/3 disks, turbulent magnetic fields become smaller than those in hot, gas pressure 
dominated disks because magnetic tension suppresses the motion with short azimuthal wave 
length. In order to show this more clearly, we decomposed magnetic fields to mean magnetic 
field B and turbulent magnetic field SB = B — B where the mean magnetic field B is 
computed by averaging the magnetic field in the region ±10 grids in zu and z directions and 
±3 grids in azimuthal direction. The bottom panels in figure 4 show the equatorial density 
and magnetic field lines depicted by the mean magnetic fields B. Figure 4e shows that as the 
disk cools, magnetic fields turn into a tightly wound spiral with well ordered mean magnetic 
fields. 

Figure 5a shows the time evolution of the magnetic pressure and Maxwell stress at 
zu = 35. The solid and dashed curves show the magnetic pressure and Maxwell stress 
normalized by the initial gas pressure averaged in the azimuthal direction and in the vertical 
direction (0 < z < 1). The magnetic pressure increases as the disk shrinks in the vertical 
direction by cooling instability. The average ratio of Maxwell stress to the total pressure 
Ob = (B^B^/An) / (P gas + -P mag ) slightly decreases but stays around «b ~ 0.1. 

Although the magnetic energy increases when low-/? disk is formed, Maxwell stress 
does not increase because the fluctuating radial magnetic field decreases. Inside the low- 
/3 disk, azimuthal component of magnetic field dominates the fluctuating radial magnetic 
field. Figure 5b shows the time evolution of the Maxwell stress normalized by the initial 
gas pressure, (B^B^/ An) / P\> averaged in azimuthal direction, in radial direction, and in 
vertical direction (0 < z < 1). Black curves show the Maxwell stress averaged in the 
inner region (4 < zu < 10), and gray curves show that averaged in the outer region (30 < 
zu < 40). The solid curves and the dashed curves show the time variation of the Maxwell 
stress computed from the mean magnetic field (B^B^/ 'Air) / and that computed from the 
fluctuating magnetic field (SB^SB^/An) / P^. In the inner region where the disk still stays in 
the hot state, the angular momentum transport is mainly due to the fluctuating component. 
In the outer region, the angular momentum transport by fluctuating field decreases with 
time and becomes comparable to that by mean magnetic fields. This figure indicates that 
the magnetic turbulence is suppressed inside the magnetically supported disk but the angular 
momentum is still transported by the Maxwell stress of mean magnetic fields. 



3.3. Increase of Luminosity 



Figure 6 shows the time variation of the X-ray luminosity computed by integrating 
the optically thin bremsstrahlung cooling term in 4 < w < 50, \z\ < 30 (solid curve) and 
in 4 < w < 50, \z\ < 1 (dashed curve) when p b = 0.29. This luminosity is normalized 
by the Eddington luminosity. X-ray luminosity increases during the growth of the thermal 
instability because the equatorial density increases 10 times while the temperature decreases 
100 times. Thus the cooling rate in \z\ < 1 becomes 10 times higher than that before the 
transition. In this stage, the effective optical depth in the \ow-f3 region is about 0.01. Thus 
the disk still stays in an optically thin state. The increase of the X-ray luminosity saturates 
because the mass accretion rate from the outer region saturates. In our simulations, the 
mass accretion rate in 30000 < t < 36000 is M ~ MEdd when pb = 0.29. Since the energy 
conversion rate r\ e from the rest mass energy to radiation is rj e < 0.06 in accretion disks 
around a Schwarzschild black hole, L/L E dd < 0.06(M/M E dd) ~ 0.06 when M ~ MecM- 
The disk luminosity shown in figure 6 is an order of magnitude smaller than this limit 
because the innermost region [w < 15) still stays in radiatively inefficient optically thin, 
hot state. The saturation level of the luminosity will increase if we start simulation with 
larger pb because the maximum accretion rate from the outer torus increases as the initial 
density of the torus increases. When the accretion rate further increases, the conversion rate 
from rest mass energy to radiation will also increase because the inner region of the disk 
undergoes a transition to radiatively efficient cool disk. In this paper, we only included the 
bremsstrahlung cooling. It should be noted that in optically thin, hot, magnetized disks, 
synchrotron Compton cooling may also contribute to increase the X-ray luminosity. 



4. Discussion 

4.1. Evolution of Optically Thin Accretion Disks 

In this paper, we have shown by direct numerical simulations without assuming a- 
viscosity that magnetically supported disk is created during the transition from a low/hard 
state to a high/soft state in black hole accretion flows. 

Figure 7 schematically shows the evolution of an optically thin accretion disk. As the 
accretion rate from the outer disk increases (figure 7a), the density of the inner accretion 
disk increases. When the density of the disk exceeds the critical density for the onset of the 
cooling instability at some radius, the disk cools and thus shrinks in the vertical direction 
(figure 7b). In such regions, magnetic pressure dominates the gas pressure because (1) gas 
pressure decreases due to cooling, while (2) the azimuthal magnetic field increases as the disk 
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shrinks in the vertical direction. The latter is due to the formation of sector structures in 
azimuthal magnetic fields. If the magnetic fields in turbulent accretion disks are purely ran- 
dom, they will be dissipated by magnetic reconnection, when the disk shrinks in the vertical 
direction. If magnetic fields have coherent structures such that the azimuthal component 
has the same sign in some regions, however, they can survive during the contraction due to 
cooling. Moreover, the field strength increases as the region shrinks conserving the toroidal 
magnetic flux. 

Since magnetic pressure supports the disk even when temperature decreases, the disk 
stops shrinking in the vertical direction and approaches a new equilibrium state. The mag- 
netically supported disk emits power-law hard X-rays because the disk still stays in optically 
thin state. Since the critical surface density for the onset of the cooling instability in the 
inner radius is larger than that in the outer radius when bremsstrahlung cooling is assumed 
(Abramowicz et al. 1995), the inward propagation of the cooling front stops and the inner- 
most region stays in the hot state. We expect that when the accretion rate from the outer 
disk further increases, the transition radius between the hot, gas pressure dominated region 
and the cool, magnetic pressure dominated region will move inward (figure 7c). 

The evolution of the innermost region of optically thin disks will depend on the cooling 
mechanism. When the synchrotron and/or Compton cooling were included, transition to the 
magnetically supported disk may start from the inner region because the cooling rate in the 
innermost region of the disk will be enhanced. 

4.2. Thermal Equilibrium Curves of low-/3 Disks 

Numerical results revealed that as the cooling instability grows, the disk approaches to a 
quasi-steady state when a low-/? disk is formed. The cooling time scale is t coo i = (P) / (Q^) ~ 
160, here ( ) represents the average of a quantity over a range of < z < 2. Figure 2b shows 
that the equatorial temperature stays T ~ 10~ 5 at least for time scale of t ~ 6000 ~ 
20tdyn- This indicates that in this state, heating balances with cooling. Since the dissipative 
heating rate of turbulent accretion disks is comparable to — {B W B V /An)m(d£l/ dm) (Hirose 
et al. 2005), the heating time scale can be estimated by thcat ~ (Pg as )/[(-B ro 5 v ,/47r)(3/2)fi], 
where f2 oc ti7~ 3 / 2 . Before the thermal instability grows, since the disk is dominated by gas 
pressure, and (B^B^/An) ~ 0.1(P gas ) (see figure 5a), t hcat = 10/[(3/2)ft] ~ 2000. After 
the cooling instability grows, the disk evolved toward a magnetic pressure dominated state 
in which (B^B^/ An) ~ 0.1(P mag ). Thus, the heating time scale at this stage is theat ~ 
(Pgas)/[0.1P ma g x (3/2)ft] ~ l/[(3/2)fi] ~ 200 at w = 35. This time scale is comparable to 
the cooling time scale. 
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Numerical results indicate that steady equilibrium solutions of optically thin, low-/? 
disks exist. Thermal equilibrium curves of accretion disks can be obtained by equating the 
heating rate Q + and radiative cooling rate Q~ d or advective cooling rate Q~ dv , assuming a- 
viscosity (Abramowicz et al. (1995); see also Kato et al. (1998)). When magnetic pressure is 
dominant, the vertically integrated pressure can be evaluated by W ~ 2B 2 H/ (8tt), where H 
is the half thickness of the disk. The vertical hydrostatic balance requires W/ (Si/) ~ Q 2 H. 
Based on the numerical results, we assume $ = BH ~ constant at a fixed radius, where 
$ is the azimuthal magnetic flux per unit radius. Thus W ~ (BH) 2 /(AtcH) ~ §1/(AnH). 
By equating this equation with W ~ EQ, 2 H 2 , we obtain H ~ [$o/(47rfi 2 )] 1/3 £~ 1/3 . So we 
obtain the heating rate 

Q + ~ ^aWn ~ ^a^ 3 {^j S 1 / 3 oc S 1 / 3 . (12) 

Since Q \ M, M \ S 1 / 3 . 

The radiative cooling rate Q; ad = Q' h p 2 T 1 ' 2 (2H) = Q'^T 1 ' 2 / \2H) is expressed as 

Q md = 9± fi 2/3 s 7/3 T l/2 K s 7/3 T l/2 

By equating Q + and Q~ d , we obtain 

When we adopt $o = {BH) 2 ~ 0.2 and a ~ 0.1 according to the results of our numerical 
simulation at t — 24100 and zu = 35, we obtain 

T ~ 0.007 x S" 4 . (15) 

Figure 8a and 8b show the thermal equilibrium curves of accretion disks. Thin curves 
show the equilibrium curves of conventional accretion disks (Abramowicz et al. (1995); see 
also Kato et al. (1998)). Thick solid curves show the optically thin low-/? branch T oc £~ 4 
depicted by using equation (15) and M oc E 1 / 3 . These curves should be connected to the 
optically thin, gas pressure dominated branch as indicated by the thick dotted curves. The 
connection between the hot gas pressure dominated branch and the magnetic pressure dom- 
inated cool branch should be studied including the advective cooling and will be presented 
in a separate paper (Oda et al. 2006 in preparation). 

Let us discuss the thermal stability of the low-/? branch. Since Q + oc E 1 / 3 (equation (12)) 
and Q- d oc Y7^T 1 / 2 (equation (13)), (<9Q+/<9T) S = and (dQ;JdT) E oc Q~ d /(2T) > 0, 

'dmQ- d \ ^ fdlnQ 

s 



8h.T ;„ > USF I (16) 



This indicates that the low-/3 branch is thermally stable (e.g., Pringle (1981)). The low-/? 
branch is also secularly stable because M oc E 1 / 3 at a fixed radius, thus the criterion for the 
secular stability (e.g., Pringle (1981)), 




Q+=Q 



rad 



> 



(17) 



is satisfied. 



Figure 8c shows the relation between E and T obtained from our numerical simulation. 
Here, T is the equatorial temperature. The symbols show the surface density and mass 
accretion rate at w = 35 averaged in azimuthal direction and integrated in \z\ < 10. The 
time range plotted is 21000 < t < 35000. An arrow indicates the time (t = 24100) when we 
included the radiative cooling. The dashed line shows the relation T oc S~ 4 expected from 
the above theory. Figure 8d shows the relation between £ and M. As the disk approaches 
the quasi-steady low-/? state, it follows the theoretically expected relation M oc E 1 / 3 . 



We would like to discuss the dependence of numerical results on the initial configuration 
of magnetic fields. Since the growth rate of MRI is larger when the initial magnetic field has 
poloidal component than the case with purely toroidal magnetic field (e.g., Hawley (2001); 
Kato et al. (2004)), the accretion from the initial torus proceeds faster in models with 
initially poloidal magnetic field. In the nonlinear stage, however, the structure and strength 
of magnetic fields do not depend significantly on the initial configuration of magnetic fields so 
long as the initial magnetic field is weak and confined inside the disk. In simulations starting 
from weak poloidal magnetic fields embedded inside the disk (e.g., Hawley (2001); Kato et 
al. (2004)), the amplification of magnetic field saturates when /? ~ 10. The magnetic field 
strength at this stage and the ratio of poloidal magnetic field to toroidal magnetic field are 
comparable to that for purely azimuthal magnetic fields. Our simulation results also show 
significant mass outflows in the coronal region, consistent with the results starting from 
weak poloidal magnetic fields. Thus, when radiative cooling is included after a magnetically 
turbulent accretion disk is formed, subsequent evolution of the disk will not be sensitive to 
the initial configuration of magnetic fields. 



4.3. Dependence on Initial Magnetic Field Configuration 
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4.4. Comparison with Observations 

Let us compare our numerical results with observations of black hole transients. Homan 
et al. (2005) pointed out that during the transition from the low/hard state to the high/soft 
state, GX 339-4 shows a sub-transition from hard intermediate state dominated by power-law 
X-ray radiation to soft intermediated state dominated by radiation from optically thick disk. 
The X-ray luminosity of the hard intermediate state is 1-10% of the Eddington luminosity 
in GX339-4 (Remillard 2005). This luminosity is much higher than the critical luminosity 
above which the optically thin, hot gas pressure dominated disk becomes thermally unstable. 
In our simulations, the critical luminosity for this transition is about 0.2% of the Eddington 
luminosity (see figure 6). The magnetic pressure dominated disk shows hard X-ray spectrum 
because it is still optically thin. Thus, numerical results suggest that the luminous disk in 
the luminous hard state stays in the thermally stable low-/? branch. 

The X-ray luminosity in our simulation results saturated at about 0.5% of the Eddington 
luminosity because the mass accretion rate from the outer torus saturated at about M ~ 
^Edd (figure 8b). We expect that when the mass accretion rate from the outer region 
continues to increase, the transition radius between the inner hot flow and the outer cool low- 
/? disk will move inward and the X-ray luminosity will continue to increase. Such simulations 
will be possible by starting simulations from a torus initially located at larger radius. 

Since the low-/? disk is supported by magnetic pressure, the transition to the optically 
thick (soft intermediate) state needs the dissipation or expulsion of magnetic fields supporting 
the disk. Recently, Fender et al. (2004) pointed out that in galactic black hole candidates, 
relativistically moving blobs are ejected during the transition from a low/hard state to a 
high/soft state. Our numerical results suggest that such ejections originate from the releases 
of magnetic energy stored in the low-/? disk. When magnetic flux inside the disk is ejected 
by such explosive events, the disk will again be supported by gas pressure and will be able 
to shrink further in the vertical direction. Extraction of magnetic flux from the low-/? disk 
enables the disk to complete the transition to the optically thick, high/soft state. 

Synchrotron and/or Compton scattering are not included in the present study. If in- 
cluded, magnetically supported disks are expected to emit enhanced radiation from the inner 
parts of the disk. It may resolve the issue of overproduction of bremsstrahlung from the outer 
parts of the low luminosity MHD flow (see Ohsuga et al. (2005) for discussion). 

In this paper, we assumed single temperature plasma. In optically thin, hot accretion 
disks, electron temperature can be lower than the ion temperature. Such two temperature 
global disk models including radial advection were studied by Nakamura et al. (1997; see 
also Manmoto et al. (1997)). To fit the observational spectra, two-temperature nature of 
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accretion flow should be considered in the simulations. Works are now in progress to extend 
our simulation code to two temperature plasmas. 
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Table 1: Units adopted in this paper. 





Numerical Unit 


Physical Unit 


length 




3.0 x 1O 6 (M/1OM ) cm 


velocity 


C 


2.9979 x 10 10 cm s" 1 


time 


u 


1O~ 4 (M/1OM ) sec 


density 


Po 


8.3 x l(r 7 (M/lOM ) g cm" 3 


temperature 


To 


1.1 x 10 13 K 


luminosity 


-^Edd 


1.25 x 1O 39 (M/1OM ) erg s" 1 
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Fig. 1. — Distribution of density (left), temperature (middle) and azimuthal magnetic field 
(right) averaged in azimuthal direction. Top panels show the initial state. Second panels 
from top show the stage before the transition (t = 24000). Third panels from top show the 
stage after cooling instability sets in (t = 32000). Bottom panels show the same stage as the 
third panels (t = 32000) for a model without including radiative cooling effect. 
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Fig. 2. — (a) Radial distribution of temperature averaged in azimuthal direction and in 
vertical direction (0 < z < 1). Curves show the temperature at t — 24000, 26000, 30000, 
and 34000 from top to bottom, (b) Time evolution of azimuthally averaged surface density, 
equatorial density, equatorial temperature, and equatorial plasma /3 at w — 35. The units 
of the horizontal axis is t — 10 4 [= 1.O(M/1OM ) sec in physical units]. 
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Fig. 4. — Distribution of equatorial density (color) and magnetic field lines projected onto 
the equatorial plane (gray curves), (top) magnetic field lines depicted by original magnetic 
fields before decomposing them into the mean field and fluctuating field, (bottom) magnetic 
field lines depicted by mean magnetic fields, (left) Before the cooling instability grows 
{t = 26000). (center) After the cooling instability sets in (t = 32000). (right) Result of a 
model without cooling (t = 32000). 
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Fig. 5. — (a) Time evolution of the azimuthally averaged magnetic pressure (solid curve) 
and the Maxwell stress (dashed curve) at w = 35 normalized by the initial gas pressure and 
averaged in < ip < 2tt and < z < 1. Dotted curve shows the ratio of the averaged Maxwell 
stress to the averaged total pressure P gas + P mag - (b) Time evolution of Maxwell stress 
computed by mean magnetic fields (solid curves) and fluctuating magnetic fields (dashed 
curves). The black curves show the Maxwell stress averaged in the inner region (4 < w < 10 
and < z < 1). The grav curves show the Maxwell stress averaged in the outer region 
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Fig. 6. — Time evolution of the bremsstrahlung luminosity integrated in 4 < w < 50, 
\z\ < 30 (solid curve) and in 4 < zu < 50, \z\ < 1 (dashed curve). 



-26- 






Fig. 7. — A schematic picture of state transitions in black hole candidates, (a) Low/hard 
state, (b) Onset of the cooling instability, (c) The transition radius between the hot disk and 
cool low-/? disk moves inward as accretion rate increases. 
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Fig. 8. — (a) Thermal equilibrium curves of accretion disks (solid curves). The thick line 
shows the optically thin, low-/? branch which is thermally stable, (b) Schematic picture of 
thermal equilibrium curves in the S — M plane. The thick line shows the optically thin, low-/? 
branch which is secularly stable, (c) The relation between surface density and equatorial 
temperature at zu — 35. The dashed line shows the theoretically obtained relation of the the 
low-/? disk T oc £~ 4 . An arrow shows the point at t — 24100 when radiative cooling term 
is included, (d) The relation of surface density and mass accretion rate at zu — 35. The 
dashed line shows the theoretically expected curve for low-/? disk M oc S 1 / 3 . 



